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Аннотация 

Введение. Широкое использование пьезоматериалов в различных отраслях стимулирует изучение их физических 
характеристик и обусловливает актуальность таких изысканий. В рассматриваемом случае модальный анализ 
позволяет определить рабочую частоту и коэффициент электромеханической связи пьезоэлементов различных 
устройств. Эти индикаторы представляют серьезный теоретический и прикладной интерес. Цель исследования — 
разработка численных методов для решения задачи определения частот резонанса в системе упругих тел. Для 
достижения цели нужны новые подходы к дискретизации задачи на основе метода конечных элементов и 
выполнение программной реализации выбранного метода на языке С# на платформе „пе. Актуальные решения 
созданы в контексте библиотеки классов комплекса АСЕГАМ-СОМРО$З. Основанные на обращении матриц 
известные методы решения обобщенной задачи на собственные значения неприменимы к матрицам большой 
размерности. Для преодоления этого ограничения в представленной научной работе реализована логика 
построения матриц масс и созданы программные интерфейсы для обмена данными о задачах на собственные 
значения с модулями пре- и постпроцессинга. 

Материалы и методы. Для реализации численных методов задействовали платформу .пеё и язык 
программирования СЯ. Валидация результатов исследования проводилась путем сравнения найденных значений 
с решениями, полученными в известных САЕ-пакетах (англ.  сотрщег-аа4е4 епошеешпе — 
компьютеризированная инженерия). Созданные подпрограммы оценивались с точки зрения производительности 
и применимости для задач большой размерности. Проводились численные эксперименты с целью валидации 
новых алгоритмов в задачах малой размерности, которые решаются известными методами в МАТЕАВ. Далее 
подход тестировали на задачах с большим числом неизвестных и с учетом распараллеливания отдельных 
операций. Чтобы избежать нахождения обратной матрицы, программно реализовали модифицированный метод 
Ланцоша. Рассмотрели форматы хранения матриц в оперативной памяти: триплеты, СВ, С$С, ЗКуПпе. Для 
решения системы линейных алгебраических уравнений (СЛАУ) задействовали итерационный симметричный 
метод ГО, адаптированный к этим форматам хранения. 

Результаты исследования. Разработаны новые расчетные модули, интегрированные в библиотеку классов 
комплекса АСЕГАМ-СОМРО$. Проведены расчеты для определения применимости различных форматов 
хранения разреженных матриц в оперативной памяти и различных методов реализации операций с разреженными 
матрицами. Графически визуализирована структура матриц жесткости, построенных для одной и той же задачи, 
но с различной перенумерацией узлов конечноэлементной сетки. Применительно к задаче теории 
электроупругости обобщены и представлены в виде таблицы данные о времени, необходимом на выполнение 
базовых операций с матрицами жесткости в различных форматах хранения. Установлено, что перенумерация 
узлов сетки дает существенный прирост производительности даже без изменения внутренней структуры матрицы 
в памяти. С учетом поставленных задач исследования названы преимущества и слабые стороны известных 
форматов хранения матриц. Так, СЗК оптимален при умножении матрицы на вектор, ЗК $ — при обращении 
матрицы. В задачах с числом неизвестных порядка 103 выигрывают в скорости итерационные методы решения 
обобщенной задачи на собственные значения. Оценивалась производительность программной реализации метода 
Ланцоша. Измерялся вклад всех операций в общее время решения. Выяснилось, что операция решения СЛАУ 
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занимает до 95 % от общего времени работы алгоритма. При решении СЛАУ симметричным методом ГО 
наибольшие вычислительные затраты нужны для умножения матрицы на вектор. Для увеличения 
производительности алгоритма прибегли к распараллеливанию с общей памятью. При использовании восьми 
потоков производительность выросла на 40-50 %. 

Обсуждение и заключение. Полученные в рамках научной работы программные модули были внедрены в пакет 
АСЕГАМ-СОМРОЗ. Оценена их производительность для модельных задач с квазирегулярными 
конечноэлементными сетками. С учетом особенностей структур матриц жесткости и масс, получаемых при 
решении обобщенной задачи на собственные значения для электроупругого тела, определены предпочтительные 
методы для их обработки. 


Ключевые слова: пьезоматериалы, метод конечных элементов, разреженные матрицы, обобщенная задача на 
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Введение. Устройства из пьезоматериалов широко используются, давно активно изучаются и 
совершенствуются. Отдельно следует отметить медицинские ультразвуковые приборы (оборудование для 
диагностики, ультразвуковые скальпели) [1-4] и мобильные генераторы энергии [5]. Исследование [6] описывает 
комбинирование фото- и пьезоэлектрического эффектов для создания действенных компактных источников 
энергии. В науке и промышленности изучаются новые материалы, рассчитанные на эксплуатацию в 
специфических условиях. В [7] рассматривается создание бессвинцового пьезоактивного состава, подходящего 
для эксплуатации при различных температурах. 

В исследовании пьезоэлементов важную роль играет этап модального анализа, позволяющий установить 
частоты резонанса и антирезонанса устройства. Эти данные: 

— необходимы для выяснения рабочей частоты устройства; 

— позволяют определить коэффициент электромеханической связи — важный индикатор эффективности 
устройства; 

— являются входной информацией в численных экспериментах для задач на вынужденные колебания. 

Цель исследования — создание численных методов решения задачи определения частот резонанса для 
системы упругих тел. Достижение заявленной цели требует решения двух задач. Первая: разработать методы 
дискретизации задачи на основе метода конечных элементов (МКЭ). Вторая: провести программную реализацию 
выбранного метода на языке С# на платформе „пе. Все известные программы учитывают контекст библиотеки 
классов комплекса АСЕГАМ-СОМРОЗ [8]. При решении обобщенной задачи на собственные значения широко 
применяются методы, основанные на обращении матриц. Однако они неприменимы к матрицам большой 
размерности. В представленной научной работе это ограничение преодолевается следующим образом: 

— дополнительно реализована логика построения матриц масс; 

— созданы программные интерфейсы для обмена данными о задачах на собственные значения с модулями пре- 
и постпроцессинга. 

Материалы и методы. В первую очередь предлагаемый подход призван решать статические задачи 
электроупругости при реализации метода осреднения [9], который задействуют для расчета эффективных 
свойств пьезокомпозитов. В связи с этим на этапе построения глобальных матриц МКЭ представлены только 
матрицы жесткости. В данном исследовании дополнительно реализовали логику построения матриц масс и 
разработали программные интерфейсы (аррИсаНоп ргоэтатите пиегасе, АР1, англ. — интерфейс прикладного 
программирования) для обмена данными о задачах на собственные значения с модулями пре- и постпроцессинга. 
Разработанные подпрограммы оценивались с точки зрения производительности и применимости для задач 
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большой размерности. Проводились численные эксперименты с целью валидации созданных алгоритмов для 
таких задач малой размерности, которые позволяют получить решение общими методами в вычислительном 
пакете МАТГАВ. Далее выполнялось тестирование на задачах с большим числом неизвестных и с учетом 
распараллеливания отдельных операций. 

Математическая модель решаемой задачи состоит из определяющих соотношений [9]: 


Р‚кб?и + азроюи -У:с0=ф, У:В=0, (1) 
с=с®.(Е+Вё)-е7.Е, О+е,О=е,: (&+6,8)+57-Е, (2) 
= = (Уи+Уи!)/2, Е =-УФ. (3) 


Здесь со — тензор напряжений; р; — плотность тела; = — тензор деформаций; и — вектор перемещений; р — 
вектор электрической индукции; Е — вектор напряженности электрического поля; /, — вектор массовых сил; 


ф — электрический потенциал; о, В.;, с, — коэффициенты демпфирования; с®, е?т, э$ — тензоры упругих 


констант, пьезомодулей и диэлектрических проницаемостей; индекс ] —Й номер тела в модели. 


Дискретизация выполняется заменой: 
и(х,г) = № (х)- 0, Фр = М(®-ФО. 


Здесь №, — матрица функций формы для поля перемещений; № — вектор функций формы для электрического 
потенциала; И (1), Ф(!) — глобальные векторы соответствующих узловых степеней свободы. 


В таком случае исходная задача (1-3) приобретает вид: 
М.-а+К-а=Е. (4) 
Здесь матрицы М и Кявляются глобальными матрицами масс и жесткости соответственно, а вектор а 
представляет собой общий вектор неизвестных: 
а= [(, Ф]. 
В задаче теории электроупругости: 


= Ми 0 К = Кии Ка 
— № о) (# -Кы/ 0 


Матрицы Мш› Кши Къ — симметричные. В случае гармонических колебаний на собственной частоте в; 


можно записать: 
а=у, 91(©и), 
обозначив через у, соответствующий собственный вектор. 
Рассмотрим свободные колебания если Г =0. В этом случае задача (4) представляется в виде: 
—0?Му +К-м, =0. (6) 

Таким образом, исходная задача сводится к обобщенной задаче на собственные значения (6). Для ненулевого 
у; неравенство (6) решается нахождением матрицы, обратной К. Однако при этом разреженная матрица 
становится заполненной, то есть метод непригоден для матриц больших размеров. Поэтому нужно использовать 
другие методы, не требующие нахождения обратной матрицы. Для решения этой задачи в данной работе 
программно реализован модифицированный метод Ланцоша [10]. Автор этой модификации — Т.С. Мартынова. 
Описание разработки в данной статье не приводится. Из используемых в методе операций наиболее затратной с 
точки зрения вычислительных ресурсов оказалось решение системы линейных алгебраических уравнений 
(СЛАУ), необходимое для выполнения спектральной трансформации. 

Матрицы М и К — разреженные, с небольшим числом ненулевых элементов. Для хранения таких матриц 
в оперативной памяти используются несколько форматов: 

— триплеты или координатный формат; 

— С5В (англ. сотргеззе4 зрагзе го\у — сжатый разреженный ряд); 

— С$С (англ. сотргеззе4 зрагзе сои — сжатый разреженный столбец); 

— формат хранилища ЗКуПпе (англ. ЗК $ — метод). 

Координатный формат предполагает хранение троек (триплетов) значений (1, }, К), представляющих собой 
координаты (1, ]) и значения (К) ненулевых элементов. СЗК иногда называют СВ$ или Йельским форматом. Он 
предполагает хранение разреженной матрицы в виде трех массивов. Рассмотрим матрицу размера № с № 
ненулевыми элементами. Опишем возможную организацию ее хранения. Все ненулевые элементы нужно 
разместить в одном массиве размера № . Позиции этих элементов в столбцах разместить в другом массиве 
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размера №7 ‚ а третий массив размера № использовать для хранения индексов первых элементов строк. 
Аналогично реализуется хранение в формате С$С. 

Формат $К$ предполагает хранение ленты матрицы переменной ширины, включающей в себя все ненулевые 
элементы. В этом случае допускается хранение нулей. Эффективность этого формата зависит от перенумерации 
строк матрицы. Методы сокращения размера ленты описаны в [11], однако требует отдельного исследования их 
применимость к матрице жесткости, получаемой при решении трехмерной задачи с использованием МКЭ. 

Для решения СЛАУ задействовали итерационный симметричный метод ГО (Зуттешле ГО Меоч, 
ЗУММГО [12]), адаптированный к перечисленным выше форматам хранения. 

Результаты исследования. В начале исследования выбрали оптимальный формат хранения для разреженных 
матриц. Координатный формат позволяет быстро добавлять и изменять элемент матрицы. Эти операции 
необходимы на этапе сборки глобальной матрицы и при учете граничных условий. Кроме того, для плохо 
обусловленных матриц, к которым относится К, часто применяют предварительное преобразование для 
нормирования. Его также удобно выполнять в координатном формате. Однако такой формат неэффективен, если 
речь идет об алгебраических операциях. 

СУК плохо приспособлен для изменения структуры матрицы: добавляя ненулевой элемент, нужно выполнять 
вставку в два массива. При этом матрица умножается на вектор очень просто и эффективно. 

ЗК$ имеет аналогичные проблемы с добавлением ненулевых элементов и сильно зависит от перенумерации 
неизвестных в задаче. Остановимся на примере квазирегулярной сетки, которая используется в пакете АСЕГАМ- 
СОМРО$З для работы с представительными объемами композитов. Ширина ленты, содержащей все ненулевые 
элементы, может быть определена заранее и зависит от числа узлов и типа конечного элемента. В общем случае 
произвольной конечноэлементной сетки сложно заранее оценить размер ленты. 

В численных экспериментах использовали четыре способа нумерации неизвестных. На рис. 1 представлена 
структура матриц жесткости, построенных для одной и той же задачи, но с различной перенумерацией узлов 
конечноэлементной сетки. 


г) 


Рис. 1. Структура матриц жесткости при различных методах нумерации узлов: а — неизвестные упорядочены по узлам; 


б — сначала упорядочены по узлам перемещения, затем потенциалы; в — узлы отсортированы по слоям КЭ-сетки, а 
неизвестные — по узлам; 6, г — узлы отсортированы по слоям КЭ-сетки, а неизвестные — как в примере 
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Итак, сетка представляла собой куб с регулярным разбиением восьмиузловыми конечными элементами. Для 
иллюстраций использовалась модельная матрица из 500 строк. Матрицы, представленные на рис. 1 а и1б, не 
подвергались дополнительной перенумерации узлов и отличаются только нумерацией степеней свободы. В Та: 

а = [Ми М.›Ф!»---›Инх>Ику >И, Фи]. 

В | бнеизвестные, отвечающие за распределение потенциала, собраны в конце вектора: 

а=[и м, И, ,,-.> Им, Иму Им ›Фрь-- Фи]. 

Неизвестные в матрицах !1в и 12 занумерованы аналогично, но узлы конечноэлементной сетки 
предварительно перенумерованы согласно их координатам путем поочередной сортировки всех узлов по каждой 
из координат. Такой прием широко используется для построения более эффективных модулей решения СЛАУ, 
так как позволяет работать с матрицей в специальном ленточном формате, удобном для распараллеливания. Для 
комплекса АСЕГАМ-СОМРО$ реализованы подобные внешние модули [13—15], однако в данной работе 
использовались только форматы хранения разреженных матриц общего вида. 

В таблице 1 сведены данные о времени, необходимом на выполнение базовых операций с матрицами в 
различных форматах. 


Таблица 1 
Время для выполнения базовых операций с матрицей жесткости в задаче теории электроупругости. 
19 652 строки 
Формат Затраченное время, мс 
Операция 

хранения Та 16 1 в 12 

С5В Конвертация из координатного формата 123 132 97 117 

С5В Умножение на вектор, 100 операций 260 260 260 260 

5К5 Конвертация из координатного формата 690 703 124 268 

5К$ Умножение на вектор, 100 операций 60 558 | 61450 7 616 22 113 


Результаты экспериментов показали, что операция преобразования из координатного формата хранения в 
компактные занимает мало времени. При этом перенумерация узлов сетки для формирования блочно-ленточной 
матрицы позволяет получить заметный прирост производительности даже без изменения внутренней структуры 
матрицы в памяти. Формат СЗВ оказался оптимальным с точки зрения эффективности операции умножения 
матрицы на вектор. При обращении матрицы формат 5К$ более эффективен, однако для задач с числом 
неизвестных порядка 103 заметно быстрее работают итерационные методы решения обобщенной задачи на 
собственные значения. 

Далее экспериментально оценили производительность программной реализации метода Ланцоша. Измерили 
вклад всех операций в общее время решения. В результате установили, что операция решения СЛАУ занимает 
до 95 % от общего времени работы алгоритма. В ходе работы алгоритма строится подпространство Крылова, и в 
зависимости от его размерности меняется число СЛАУ, которые необходимо решать. Отметим, что размерность 
подпространства Крылова выбиралась на основе эвристик относительно числа искомых собственных значений. 
При этом СЛАУ отличаются только правой частью, благодаря чему требования к выделяемой памяти остаются 
невысокими. Среди базовых операций, применяемых в ходе решения СЛАУ методом ЗУММГО, наибольшие 
вычислительные затраты нужны для умножения матрицы на вектор. 

Для увеличения производительности алгоритма реализовано простейшее распараллеливание с общей 
памятью. Для формата СВ$ выделялись блоки строк. Они передавались в отдельные потоки, которые вычисляли 
соответствующие компоненты результирующего вектора. Прирост производительности составил 40-50 % при 
использовании 8 потоков. При этом для матриц порядка 103 элементов прирост был около 40 %, а для матриц 
порядка 10% — около 50 %. 

Обсуждение и заключение. В рамках данной работы реализован метод решения обобщенной задачи на 
собственные значения для матриц, получаемых при моделировании электроупругих тел. Созданы программные 
модули на языке С# для построения матриц масс методом конечных элементов и выполнения вспомогательных 
операций в рамках метода Ланцоша (работа с векторами подпространства Крылова, переортогонализация, 
нахождение собственных векторов). Вычислительная сложность обусловлена в основном операциями 
умножения разреженных матриц на вектор. В связи с этим проводились численные эксперименты по 
определению оптимальных форматов хранения матриц, оптимальной структуры матрицы, получаемой в 
результате перенумерации узлов КЭ-сетки и степеней свободы в узлах. Разработана версия итерационного 
алгоритма ЗУММГО с использованием параллельных вычислений. Итоговая схема работы включает три пункта. 


Оганесян П.А. и др. Реализация базовых операций для разреженных матриц 


Во-первых, строятся глобальные матрицы в координатном формате с алгоритмом перенумерации (рис. 1 в). Во- 
вторых, данные преобразуются в формат СЁ$. В-третьих, данные обрабатываются методом Ланцоша, который 
включает метод ЗУММГО для решения СЛАУ. Результаты работы включили в программный пакет АСЕГАМ- 
СОМРОЗ. 
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